{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy import constants as const\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"计算总的UV流量，sumUV为21cm连续谱的总数（范围：l，b=51.892，0.562 半径14.8角分），单位为K。\n",
    "1pixel = 18arcsec。根据换算关系5.5my=1K（https://www.ucalgary.ca/ras/VGPSintroduction）\n",
    "\"\"\"\n",
    "def getSv(sumUV):\n",
    "    beam=60*u.arcsec\n",
    "    pix = 18*u.arcsec\n",
    "    total_flu = sumUV*5.5*u.mJy/(beam/pix)**2\n",
    "    return total_flu.to(u.Jy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"根据柱密度计算质量，1pixel = 3arcsec\n",
    "Nh2 = 3.3691056e+27\n",
    "npix = 426138\n",
    "Nh2bg = 5.4629e+21（SumNh2Bg=2.3279492802e+27） 可知背景的噪声是非常大的，对结果有很大影响\n",
    "\"\"\"\n",
    "def totalMass(sumNh2):\n",
    "    mu=2.8\n",
    "    pixSize=((3/3600*np.pi/180*9.63*1000*u.pc)**2).cgs\n",
    "    totalMas=pixSize*sumNh2*u.cm**-2*mu*const.m_p\n",
    "    return totalMas.to(u.M_sun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"计算平均密度\n",
    "\"\"\"\n",
    "def nmean(Mass,R):\n",
    "    mu=2.8\n",
    "    n0 = (Mass*u.M_sun/(4/3*np.pi*R**3*u.pc**3)/mu/const.m_p).cgs\n",
    "    return n0   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"返回银心距离\n",
    "\"\"\"\n",
    "def getRg(ll,bb,vlsr):\n",
    "    Rg = 8.5*np.sin(ll*np.pi/180)*(220/(vlsr+np.sin(ll*np.pi/180)*220))\n",
    "    return Rg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rg = getRg(51.892,0.562,4.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.29865402660295"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"计算电子温度,Te_error = 350+46*Rg\n",
    "\"\"\"\n",
    "def getTe(Rg):\n",
    "    Te = 5780+287*Rg\n",
    "    return Te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8161.713705635047"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getTe(Rg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "731.7380852237357"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "350+46*Rg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "R=41.5\n",
    "Sv=getSv(77769.678)\n",
    "mass = totalMass(3.3691056e+27)\n",
    "n0 = nmean(mass.value,R)\n",
    "Rg = getRg(51.892,0.562,4.2)\n",
    "Te = getTe(Rg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sv= 38.49599060999999 Jy &mass= 1482212.3831570582 solMass \n",
      " n0= 71.54445771984301 1 / cm3 &Te= 8161.713705635047\n"
     ]
    }
   ],
   "source": [
    "print(\"Sv=\",Sv,\"&mass=\",mass,\"\\n\",\"n0=\",n0,\"&Te=\",Te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8161.713705635047"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tef = Te+"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv=1.4\n",
    "D=9.63\n",
    "Nly = 4.76e48*mv**0.1*Te**(-0.45)*(Sv)*(D)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$3.0520119 \\times 10^{50} \\; \\mathrm{Jy}$"
      ],
      "text/plain": [
       "<Quantity 3.05201188e+50 Jy>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Nly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Nly=Nly.value\n",
    "aB=(10000/Te)**-0.75*2.6e-13\n",
    "Rs=(3*Nly/((4*np.pi)*n0**2*aB))**(1/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs=Rs.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "R=((41.5*u.pc).cgs).value\n",
    "Cs=((10*u.km/u.s).cgs).value\n",
    "#动力学年龄\n",
    "tr=(4/7*Rs/Cs*((R/Rs)**(7/4)-1)*u.s).to(u.year)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$4831453.9 \\; \\mathrm{yr}$"
      ],
      "text/plain": [
       "<Quantity 4831453.92766677 yr>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"碎裂时标\n",
    "\"\"\"\n",
    "tf=1.56*(100/1000)**(-5/11)*(Nly/10e49)**(-1/11)*1**(7/11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$4.0143641 \\; \\mathrm{\\frac{1}{Jy^{1/11}}}$"
      ],
      "text/plain": [
       "<Quantity 4.0143641 1 / Jy(1/11)>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
